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Tensor network states and specifically matrix-product states have proven to be a powerful tool for simulating 
ground states of strongly correlated spin models. Recently, they have also been applied to interacting fermionic 
problems, specifically in the context of quantum chemistry. A new freedom arising in such non-local fermionic 
systems is the choice of orbitals, it being far from clear what choice of fermionic orbitals to make. In this work, 
we propose a way to overcome this challenge. We suggest a method intertwining the optimisation over matrix 
product states with suitable fermionic Gaussian mode transformations. The described algorithm generalises 
basis changes in the spirit of the Hartree-Fock method to matrix-product states, and provides a black box tool 
for basis optimisation in tensor network methods. 


Capturing strongly correlated quantum systems is one of 
the major challenges of modern theoretical and computational 
physics. Recent years have seen a surge of interest in the 
development of potent numerical methods based on tensor 
networks to approximate ground states of interacting lattice 
models GH2 , building upon the success of the density-matrix 
renormalisation group (DMRG) (T). It has become clear that 
such ideas are also applicable to fermionic systems EHEI. 
and even to systems of quantum chemistry UBED, lacking 
the locality present in lattice models in condensed-matter sys¬ 
tems. Such tools allow in principle to approximate the full 
configuration interaction solution to good accuracy with rea¬ 
sonable effort, going in instances beyond conventional ap¬ 
proaches to quantum chemistry, such as coupled cluster l20l . 
configuration interaction or density-functional theory EDG2, 
as convincingly shown by first implementations of DMRG al¬ 
gorithms in quantum chemistry (QC-DMRG) fTTHl5l . 

Yet, there is a new obstacle to be overcome: Tensor net¬ 
work methods have originally been tailored to capture local 
interactions, and consequently ground states exhibiting short- 
range correlations and entanglement area laws JT). Systems 
in quantum chemistry pose new challenges due to the inherent 
long-ranged interactions, which are present no matter in what 
basis the systems are expressed. New questions hence arise 
concerning the optimal topology and physical (orbital) basis 
used to construct the tensor network state dGHElMl. 

In this work, we propose a novel approach towards mak¬ 
ing use of tensor network methods in quantum chemistry, 
by suggesting an adaptive scheme of updating basis transfor¬ 
mations “on the fly” in conjunction with tensor network up¬ 
dates. In this way, we bring together advantages of matrix 
product states - which can capture strongly correlated states, 
but are tailored to short-ranged correlations and low entangle¬ 
ment - and fermionic Gaussian mode transformations - for 
which entanglement is no obstacle, but non-Gaussian corre¬ 
lations are. We hence go significantly beyond previous ap¬ 
proaches towards optimising fermionic bases in tensor net¬ 
work approaches to quantum chemistry. Previous DMRG im¬ 
plementations in quantum chemistry allowing for an optimisa¬ 


tion of the physical basis restrict the mode transformations to 
permutations and separate the optimisation over the basis and 
state such that multiple DMRG runs are necessary |26). As 
a first attempt basis optimisations using a few transformations 
have been implemented for tree tensor networks, however, this 
has been found to be unstable ED. Mixing fermionic orbitals 
from an active space - the space considered here - with fur¬ 
ther ones from an additional external space has also been stud¬ 
ied H27I429I . In these approaches orbital transformations are 
carried out again between different DMRG runs. In contrast 
to this we perform the mode transformations within the active 
space in parallel to the state optimisation and directly optimise 
the entanglement structure of the tensor network. 

We focus on matrix-product states, but explain in what way 
the idea is generally applicable. We also discuss the role of 
symmetries and the geometry of the problem at hand. The ba¬ 
sis optimisation is incorporated into the standard two-site QC- 
DMRG and can be added to existing implementation without 
increasing the computational costs of the DMRG. The result¬ 
ing scheme can be used in parallel to a ground state search 
or as a pre-processing step in which the physical basis is op¬ 
timised in a first phase restricting the bond dimension of the 
MPS used to medium values and calculating the final ground 
state in the optimised basis with higher accuracy. 

System. In this work, we are concerned with strongly cor¬ 
related interacting fermionic models with a finite number of 
relevant modes as they appear in the quantum chemistry con¬ 
text. In second quantised form the Hamitonian takes the form 

np np 

H = ^2 Ujcjcj + ^2 Vi,j,k,ic\c]cic k , ( 1 ) 

i,j =1 i,j,k,l =1 

where Cj is a fermionic annihilation operator associated to the 
mode labeled j satisfying the canonical anti-commutation re¬ 
lations {ci,cj} = 0 and {c\,Cj} = e> t j and the coupling t 
and v are such that H is Hermitian. p denotes the number of 
different fermion species present for each of the n orbitals, 
e.g. spin up and down electrons. The one particle modes 
form the basis of single particle Hilbert space 'H rip . Any 
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fermionic state will be an element of the fermionic Fock space 
F = A A di np , where A denotes the exterior product 

and A °'H np = C, with a basis formed of all Slater deter¬ 
minants |cc), where x £ {0, l} np , of the initial single parti¬ 
cle modes. We refer to this basis as the physical basis. The 
Jordan Wigner transformation establishes an isomorphism be¬ 
tween F and the Hilbert space of n qudits "H®" = C d " with 
p = log 2 d. By choosing any ordering of the orbitals such 
systems can be viewed as one-dimensional lattices of n sites 
with long-range interactions. 

MPS and general idea. For a one-dimensional quan¬ 
tum lattice with n sites, where each site is described by 
a d— dimensional Hilbert space PLd, a matrix product state 
(MPS) vector takes the general form 

d 

M = E A ll) ■■■*$] (2) 

ai,...,ot n =l 

where A£ (;Dm-ixD m and {| a )} form a basis of T-L d and 
Do = 1 = D„. If the bond dimension 1) £ N is allowed 
to vary arbitrarily over different sites, every quantum state 
of the lattice can be written as in Eq. 0 130). Restricting 
the maximal bond dimension along the chain to a fixed value 
D max creates the sub-manifold Ad £> max of the full state space. 
Approximations of the ground state of a given Hamiltonian 
within this sub-manifold can be found using the density ma¬ 
trix renormalisation group algorithm (DMRG) which, as an 
alternating least square method, optimises the entries of the 
MPS tensors (A[ m ]) m iteratively E5ll3Bl33l . 

The freedom one has in this construction is a redefinition of 
the fermionic modes by a linear transformation. Linear trans¬ 
formations of a set of fermionic annihilation operators { c ,} 
to a new set {di} satisfying the canonical anti-commutation 
relations are captured by Cj = Y^j=i Uijdj, with a unitary 
mode-transformation U £ U(np). This change of the single 
particle modes induces a transformation of the physical ba¬ 
sis of F. Under this change of basis a fermionic state vector 
|^>(1)) transforms to \if(U)) = with the Gaus¬ 

sian unitary transformation G(U) = exp []TA . (In U‘ { )i,jc\cj\ 
acting in Fock space. The transformation on C d induced 
from the Jordan Wigner transformation is given by g(U) = 
©fc=o A fe U* where A° & = 1. 

We now turn to describing ground states of fermionic 
Hamiltonians with MPS expressed in a given basis, where 
the approximatability of the states strongly depends on the 
choice of basis mm. Specifically, denoting the Hamil¬ 
tonian written in terms of the transformed modes by H (U ) = 
G(Uy HG(U), we are interested in the solutions of 

(U 0 P t, IV'opt)) = argmin {^\H(U)\ip). (3) 

ueu{np),\^)eM Dam 

Note that the Hartree-Fock method is readily included in Eq. 
0 by the case D max = 1, when \if) is restricted to be a 
Slater determinant. Identifying the optimal or close-to opti¬ 
mal basis for a general Hamiltonian and D max in the sense 



FIG. 1. Illustration of the general ansatz-class of an MPS with vary¬ 
ing physical basis, where g(U) is a Gaussian transformations defined 
by a mode transformation U £ U ( np ) as described in the main text. 

of Eq. 0 would provide a deeper understanding of the en¬ 
tanglement structure of ground states appearing in quantum 
chemistry, but since this is a non-convex problem, approxi¬ 
mate solutions are accessible only. Here, we take an approach 
that iteratively finds close to optimal solutions numerically by 
optimising over the ansatz-class depicted in Fig. |T] 

Compositions of local mode transformations. In order to 
calculate approximations to the solutions of Eq. 0 and avoid¬ 
ing stability and performance issues of a direct global opti¬ 
misation, we perform successive local mode transformations 
in parallel to a two-site QC-DMRG and use a few additional 
global reorderings of the orbitals as in Refs. [24] (26) to leave 
local minima during the optimisation-process. Given a state 
vector a site-index m £ [n — 1] and a cost function f m 
which will be discussed below we solve 

C/]°] = arg min f m (\ip(l pm © U © l pn -pm- 2 p ))), (4) 

uev 

with Ifc denoting the k —dimensional identity matrix and V C 
U(2p) needs to be chosen depending on the symmetries of 
the system. The global basis change is then composed of lo¬ 
cal unitaries which are solutions of 0 for different m and act 
non-trivially on overlapping areas of the lattice and interme¬ 
diate global reorderings of the lattice-sites. 

The cost function is chosen according to the following 
paradigm. The bond dimension needed for a bipartition of the 
system to approximate a state up to a predefined accuracy can 
be upper bounded using the Renyi entropies S a (p re d) of the 
reduced state for a < 1 l34l . where S a (p) = logtrp Q /(l — 
a). We therefore iteratively minimise the ,S'i entropy over 
the chosen bipartition by using the cost function /m' > (|V’)) = 
||£™||i where £™ denotes the Schmidt spectrum of | if)) for 
a bipartiting cut between sites to and to + 1. With increas¬ 
ing dimension of V, which growth with the number of species 
per orbital p, and bond dimension D max the optimisation of 
fm' 1 becomes slow. Efficiency can be gained by minimis¬ 
ing /A 1 '* = — 1|£™||4 of which we can calculate the gradient 

vUij fnV(\ip{lpm © u ® lpn-pm- 2 p))) analytically and effi¬ 
ciently in the bond dimension as shown in the appendix. The 
optimisation of S 2 will not lead to certified bounds on the re¬ 
quired bond dimension, but will favour stronger decays in the 
Schmidt spectrum similar as the minimisation of S'?. 

The results presented here have been obtained by optimis¬ 
ing the one norm of the Schmidt spectrum, fm\ The optimi¬ 
sation of fmf can be applied in if V has a higher dimension; 
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which appear for p > 1 if the system lacks symmetries. Both 
the choice of the cost function and symmetries influence the 
choice of V as argued in the following. 

Optimisation set. In the presence of symmetries, choos¬ 
ing a physical basis which can be labeled by good quantum 
numbers decouples different symmetry sectors in the coef¬ 
ficient tensors of the Hamitonian and the MPS and allows 
for more efficient computations. Only mode transformations 
which commute with the generators of the symmetry transfor¬ 
mations will preserve the structure imposed by the symmetry. 

In general QC-DMRG algorithms only exploit a subgroup 
of the full symmetry group of a specific Hamiltonian, such as 
conservation of the number of particles, spin reflection sym¬ 
metries, Abelian point group symmetries or a 517(2) spin ro¬ 
tation symmetry l[Tl fl6l[35T4381 (see also Ref. Ii25l ). Here, we 
consider for the states the case of particle number conserva¬ 
tion of each species, which is an Abelian symmetry and allows 
for an easy implementation of symmetric MPS j39l and want 
the local mode transformations to respect the SU (2) symme¬ 
try of the considered systems. The admissible transformations 
in this case are of the form U = U® p with U n £ U (n) acting 
on one species of fermions. 

The cost functions f m chosen above depend only upon 
the Schmidt-spectrum of the state for a cut between sites m 
and to + 1 and are therefore insensitive to mode transfor¬ 
mations of the form U m © (7 m +i with U q £ U{p) acting 
only on the modes associated to the lattice site q. To obtain 
a non-redundant parametrisation of the unitaries used in the 
optimisation in Eq. 0 we restrict it to the set of left cosets 
U{2p)/U{p) x U(p) which is isomorphic to the Grassmann 
manifold G(2p,p). Efficient implementations of optimisa¬ 
tion algorithms such as the conjugate gradient method within 
Grassmann manifolds using 2 p 2 parameters are described in 
Refs. l40l BTI If we restrict ourselves to mode transforma¬ 
tions which preserve the ,S77(2)-symmetry. the relevant mode 
transformations are parametrised by G(2,1), leaving 2 opti¬ 
misation parameters in each step. Focussing on this case here 
with medium values for £> max , we can obtain U ] °[ of fin 1 by 
using gradient free schemes such as the Nelder-Mead method, 
due to the small number of parameters. 

Algorithm. Combined with an approximation of the ground 
state of a given Hamiltonian, local mode transformations nat¬ 
urally extend a two-site DMRG. A single two-site DMRG step 
results in a blocked tensor A[ m>m+1 ] £ (£d 2 x-D m _ix.D m+ i. 
In the generic case, restoring the MPS format in Eq. 0 by 
decomposing the blocked tensor into local tensors £ 

C dxD m _ lX D; and A[ m+1] £ c d *D' m xD m+1 will lead to 

D' m > D max such that the found state needs to be projected 
into A4 £> max by discarding the D max — D' m smallest values of 
the resulting Schmidt spectrum E™. The projection yields an 
truncation error = Yii > where m £ E™ are the dis¬ 
carded singular values. If we allow for a local mode transfor¬ 
mation before the blocked tensor is decomposed, the trunca¬ 
tion error can be reduced. 

Using the gauge-invariance of MPS, we bring the MPS 


TABLE I. Two site DMRG with adaptive mode transformations. 
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2 

3 

4 
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iterate over neighbouring sites m £ [n — 1]: 

get blocked tensor A ^ m+1 , (e.g. from two-site DMRG) 
calculate (local) minimum U \ p j of / m (|t/>(l © U © 1)}) 

if /mMl ® U‘° p j © 1)} < /m(|V>(l)»: 

transform | ip) with U by updating 


_ sr^d /rrloc\ jc/,/3 7 

[m,m-\- 1] 2-^oc',/3'= 1 opty(a,/3),(a-',/3 , )^[ mjm+ i] 

transform relevant operators with t/' 


and 


calculate A [m j, A [m+1] by decomposing A [mim+1) with 
truncation and update MPS with new tensors 


to a mixed normalised form, i.e. matrices of sites q < m 
are left-normalised whereas matrices associated to sites q > 
in + 1 are right-normalised j2l |42| . We can then calculate 
E™ from the blocked tensor A[ miJn+1 ]. We optimise the ba¬ 
sis by solving Eq. 0 while keeping expectation values of 
the state, such as the energy, constant by using (ip\H\il)) = 

with U = l pm © © v_ pm _ 2p . 

As the mode transformation acts non-trivially only on sites to, 
to + 1 the transformed state vector | i/j) can be represented by 

A[k](U) = A[ k ]( 1), k £[n]\{m,m + 1}, (5) 

A [im + ^U)= T, 9(U) {aMa £^)A^ +1] (l). ( 6 ) 

a',/3' 

Operators such as the Hamiltonian can be trans¬ 
formed efficiently using their second quan¬ 
tised representation. For an operator O = 

— 1 ‘ C l S C jt • ‘ ■ C jl With 

o = o(l), the coefficients transform under a mode- 

transformation according to o(U) = (W)® s o(l)U® r . As 
most operators of interest, e.g., the Hamiltonian, contain 
terms with small s and r those transformations can be 
implemented efficiently; with cost scaling as 0((np) s+r_1 ) 
for local transformations. 

Standard QC-DMRG algorithms use complementary oper¬ 
ators um m] Q5i in order to reduce the computational cost 
of each DMRG step. In the appendix we show that comple¬ 
mentary operators transform as general operators under local 
mode transformations and argue that local mode transforma¬ 
tion can be found and applied in a time not exceeding the com¬ 
putational cost of a single DMRG step. This allows us to keep 
the structure and the computational complexity of the two-site 
DMRG algorithm and perform basis optimisations essentially 
for free with the algorithm in Table |T] 

Numerical results. We use a QC-DMRG algorithm which 
uses the dynamical block state selection approach m and 
configuration interaction based dynamically extended active 
space m procedure to accelerate the convergence and adapt 
the basis of the physical space by the algorithm described 
above. As a test-system, we have chosen the electron- 
configuration of a Beryllium ring built from 6 Be atoms. This 
system has recently been investigated l24l and a strong depen¬ 
dence of the convergence of the DMRG from the initial basis 
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FIG. 2. Numerical results for the Be ring of 6 Be atoms with in¬ 
teratomic distance of 3.3 A. All calculations have been performed 
with a 17(1) x (7(1) symmetric open boundary MPS and local mode 
transformations which keep the 517(2) symmetry of the Hamilto¬ 
nian. Both diagrams show results of the described optimisation for 
the physical basis. In the left panel we show the bond dimension 
needed for a bounded truncation error e t rc < 10 -6 and Dmin = 64 
when starting in the HF basis. The dark blue dashed line corre¬ 
sponds to a calculation in the HF basis, the blue dotted and light 
blue line correspond to the first and the tenth iteration of the calcu¬ 
lation with basis optimisation. The right panel compares the relative 
error in energy ((ip\H\ip) — Eq)/Eq obtained by calculations with 
D mdX = 256, where the reference value for the ground state energy 
Eo has been obtained from a calculation with D mdx = 2048 in the lo¬ 
calised basis. The dark blue dashed and red dashed-dotted line show 
the results for a calculation in the HF and localised basis respectively. 
In light blue and orange we plot the relative error of the 15th and 10th 
iteration of the calculation with basis optimisation starting in the HF 
and localised basis respectively. 


was observed. We investigate the molecule in a stretched ge¬ 
ometry with an interatomic distance of 3.3 A. As initial bases 
we use the Hartree-Fock (HF) basis of the system and a lo¬ 
calised basis derived from the HF orbitals by a Foster-Boys 
localisation ED- Such localised orbitals are widely used in 
QC-DMRG calculations and are known to yield a better con¬ 
vergence for the Be ring lf24l . 

Starting from the according initial basis we iteratively ap¬ 
ply the following scheme: we run the standard QC-DMRG 
for 2 sweeps, perform 8 additional sweeps together with the 
local mode transformation as described in Tableland reorder 
the basis according to its mutual information patterns j26| . 
Hereby we either fix the truncation error e t rc made in each 
step or set a hard-cut on D max . The results of our calculations 
are show in Fig. [2] In the left panel of Fig. [2] we show how 
the bond dimension behaves for a ground state search with a 
bounded truncation error e trc < 10 fi and /7 lmn = 64 for a 
calculation in the HF basis and optimised bases obtained by 
the above scheme starting from the HF orbitals. 

It is key to the method proposed that the optimisation of 
the basis leads to a significant decrease of needed resources 
already in the first iteration, where after the tenth iterations 
of basis optimisation the needed bond dimension is more than 


one order of magnitude smaller than in the unoptimised or¬ 
bitals. For realistic applications of the scheme intermedi¬ 
ate high bond dimensions during the calculation need to be 
avoided. The right panel of Fig. [2] shows the relative er¬ 
ror in energy reached when performing a calculation with 
D max = 256 starting in the HF and localised basis. As noted 
before the localised basis allows for a more efficient approx¬ 
imation of the ground state than the HF basis. The basis op¬ 
timisation allows us to further significantly optimise both the 
HF and the localised basis. Starting from the HF orbitals we 
obtain a basis for which the relative error in energy drops by 
one order of magnitude from 1.52 x 10 -4 to 1.2 x 10 -5 . Be¬ 
ginning at the localised basis allows us to reduce the relative 
error in energy from 3.7 x 10 -6 to 8.3 x 10 -7 . Note that 
the energy in the optimised basis starting from the HF orbitals 
is slightly worse than the energy obtained in the localised or¬ 
bitals, reflecting the fact that finding the optimal basis is a hard 
global optimisation problem. 

We have repeated similar calculations for different config¬ 
urations of the Be ring; at the equilibrium configuration and 
close to the avoided crossing. In each case we have been able 
to find a physical basis allowing for a more efficient approx¬ 
imation of the ground state. This illustrates that the above 
scheme can significantly and efficiently optimise a given ini¬ 
tial basis. As the local mode transformation can be added with 
no increase of the computational cost to an existing two-site 
DMRG and typically yield already in the first iteration of the 
basis optimisation a significant improvement of the basis our 
scheme extends the standard QC two-site DMRG. 

Conclusion and perspectives. In this work, we have pre¬ 
sented a scheme that adapts the physical basis an MPS is for¬ 
mulated in by applying Gaussian transformations. Incorporat¬ 
ing local Gaussian transformations into the two-site DMRG 
algorithm allows us to optimise both the basis and the MPS 
iteratively. The resulting algorithm successfully optimises the 
physical basis such that distinctly better approximations of the 
ground state of a given system by an MPS can be identified. 

It should be manifest from the description of the method 
that the same idea is equally applicable to other tensor net¬ 
works, due to the locality of the transformations. In particu¬ 
lar, tree-tensor network approaches IT9ll44il45l can readily be 
combined with the methods laid out here. Similarly, they are 
expected to be helpful for 2-d lattice systems ||3lf4l[3~2|. In ad¬ 
dition the above scheme can be directly combined with recent 
developments for the time-evolution of MPS li46l in order to 
obtain a time-evolution with variational physical basis. 

Our general strategy - of combining tensor networks with 
fermionic transformations - complements the recent interest¬ 
ing approach of Ref. B2. which is similar in mindset, but 
where these two components are put together in the oppo¬ 
site order. There, a matrix-product operator is applied onto 
a free fermionic wave function. In contrast to that approach, 
we here retain efficient contractibility, however. The approach 
taken in this work can also be seen as a variational princi¬ 
ple that allows to find the optimal fermionic tensor network 
in Ref. 1481 . where a fixed fermionic basis change is being 
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made use of. Widening the scope, these tools seem also help¬ 
ful in related approaches making use of a big data machinery 
to capture strongly correlated quantum systems. For example, 
compressed sensing ideas can help finding localised Wannier 
functions |49 50|. which in turn can be made use of in den¬ 
sity functional theory ll2Tl [22l . In conjunction with the tools 
developed here, a combined approach close to optimally rep¬ 
resenting fermionic correlated states seems within reach. 
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APPENDIX 

Gradient and geometry of the optimisation problem 

Implementing the optimisation methods in Grassmann 
manifolds as described in Refs. mm we parametrise the 
Grassmann manifold G(a, b) by a x b isometries X £ C axb 
which form the Stiefel manifold V(a. h). In order to imple¬ 
ment a the conjugated gradient search for identifying a (lo¬ 
cal) minimum of / : G{a , b) —> R. in G(a, b) the derivatives 
df(X)/dReXij + idf(X)/dImXij are needed. As the re¬ 
sulting mode transformation acts on two sites only, we are 
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with 


r r _ 

H '° = da ’ 

Zi = (X — P)(l| — A^P) -1 , 
Z 2 = (la - - pt). 


( 12 ) 

(13) 

(14) 


Preservation of the DMRG structure 


FIG. 3. Tensor network representing ||S^ lm 0 [ / 0 1 _ m _ 2 )||t with 
U G U (2 p) if j ip) is represented by a MPS in mixed normalised form 
as discussed in the main text with being the blocked tensor 

of sites m and m + 1. Legs corresponding to indices of site m are 
indicated in blue, to indices of site m + 1 in red. 


interested in the case b = a/2. 

Given an isometry X G C axa / 2 , i.e. X^X = 1° we can 
construct a unitary U(X) G U(n ) with the first a/2 columns 
being equal to the columns of A' by 

U{X) = 1„ - (X - P)( 1- - Atp)- 1 (At - pt), (7) 

with P G C“X“/ 2 and P t ) = (/j which corresponds to a 
generalised Householder reflection of the subspace spanned 
by the columns of X. Note that if l* — X'P turns out to 
be singular, we can always transform the columns of X using 
a random a/2 x a/2 unitary, e.g. with random, by 

which we chose a new representative in V(a, a/2) of the same 
element in G(a, a/2 ) 

Given a general invertible matrix Y, the elements of which 
depend on a real parameter t, we can evaluate the derivative 
of the inverse matrix to 

^Y®- 1 =Y(t)- 1 [±Y(t)]Y(t)~ 1 . (8) 

From Fig. [3] we can read off 

= tr([ 5 (P)t 0 g(U)']M\g(U) 0 g(U)]N), (9) 

where M(A[ m rn+1 1 ) corresponds to the inner ring of tensors 
in Fig. [3] which depends on the coupled tensor ,4r m m+1 ] and 
N orders the outer legs as shown. We can then evaluate the 
derivative of this cost function with respect to the parameters 

Re(Xi } j) and Im(Xij) using 


dg{U)i,j 

dU h 


( _ ir w+P,(i) detf/ t| A{ . }iA0 . } 


GO) 


if |/| = \J\,i G I and j G J and 0 otherwise, where U'\i j = 
(Ui j)izi,j£j, I,JC [np\ and px(x) denotes the number of 
elements of X smaller x and 


dU(X(a)) 

da 


-X’Z 2 - Z x X* + Z 1 X^'PZ 2 , 


( 11 ) 


To avoid redundant calculations, DMRG implementa¬ 
tions use complementary operators which are expanded and 
reloaded during the sweeps mmm. We denote by L m the 
set of all modes which are associated to sites q with q < m and 
by R m the set of modes that belong to sites q with q > m +1. 
In addition we define the abbreviations 

d 

\ L -)= Y, (4lf 1 )...4l“™_- 1 1 ] ) a |a 1 )0...0|a m _ 1 ) 

(15) 

and 


I Ra)= Y (^]...A“l) 0 |a ro+1 )®...®|a„). 

Ct-m +lv-5 Q: n, = l 

(16) 

The four different types of complementary operators, 

P\j a '\ Q T C k ’ a ’\ sY a ’ b with I m = L m , R m , a, b G 

[D m _ i] and i,j, k G [np]\{I m } are defined as follows 


P i Lm idib , 

i,j ~ t r [m-l] 


(\L a )(Lb\ 0 1 ) Y v i,j,k,i4 c j c i c k 


k,l£Lr, 


(17) 


Qi,k’ a ’ b = tT [m- 1 ] (\L a )(L b \ 0 1 ) Y V i,j,k,l C Wj ClCk 


T~)I J m ,CL,b , 

R, = trr m _i 


oLm i a ,b _ 


(18) 

(\L a )(L b \ 0 1) Y v i,j,k,l4 C j ClCk 

j,k,leL m 

(19) 

(\L a )(L b \®t) Y kAcj], (20) 

j£L m 


where the partial traces are evaluated according to the Jordan 
Wigner representation, plus their corresponding counterparts 
on sites q with q > m + 1 which are defined analogously by 
replacing L m by R m , (\L a )(L b \ 0 1) by (1 0 |i? a )(i? b |) the 
partial trace over sites [m — 1 ] by the partial trace over sites 
[n]\[m + 1], During a DMRG run the complementary oper¬ 
ators for different m are saved and reused in later steps. In 
order to update the complementary operators efficiently with¬ 
out loading and saving many operators from and to the disk, 
we evaluated some of the basis changes in a lazy fashion. Dur¬ 
ing a right sweep, the complementary operators l >R,n , Q Rm , 
R r ™ and S R ™ are loaded while a mode transformation trans¬ 
forming the orbitals 1,..., m — 1 was accumulated (resulting 
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from the previous 2 m — 3 steps). We obtain the complemen¬ 
tary operators within the updated basis by transforming them 
similar to general operators as discussed in the main text by 


{U)= Y, 

: a>6 (l) 

i' ,j'£[np]\Im 


= [{U\®u\)P L ^ b (l)]. l p, 

(21) 

(U) = £/jQ J -“> 6 ( l)E/>, 

(22) 

(U) = UjR Im ’ a ’ b ( 1 ), 

(23) 

(U) = u}S Im ’ a ’ b {l), 

(24) 


with I L, II for a left and right sweep respectively and 
Ul and Ur the corresponding accumulated transformations. 
Note that, for a right sweep, the operators P im , Q Lm , R Lm 
and S Lm can be formed with the rotated coefficients of the 
Hamiltonian and need no further rotation. 


Cost-analysis for the presented scheme 

The run time for a two-site DMRG per DMRG-step in¬ 
cluding an update of the complementary operators when step¬ 
ping from site to to m ± 1 during a sweep as 0(n 2 D 3 2 2p + 
n 2 D 2 2 3p + n 3 D 2 p 2 ) lfl4ll . The steps of the algorithm pre¬ 
sented in Table 1 in the main text scale in addition to that as 
follows. Note that p is typically 2 (spin up and down elec¬ 
trons) for a QC-DMRG - but we will keep track of the de¬ 
pendence of the cost on p for completeness. The rotation 
of the blocked matrix and computation of f„, comes at a 
cost of (){I)‘ 2 2 Ap + D 3 2 3p ). The transformation of the con¬ 
jugated operators are performed with a cost of 0{n 3 p 3 D 2 ) 
for general transformations, where the cost can be lowered to 
0(n 3 p 2 D 2 ), once a U(l) xp or more general symmetry is en¬ 
forced on the transformations. Due to the locality of the mode 
transformation, we can transform coefficients of the Hamilto¬ 
nian in a time scaling as 0(n 3 p 3 ). The computational cost 
induced from the local mode transformations which optimise 
f A n iteratively are therefore given by 0(D 2 2 4p + D 3 2 3p + 
n 3 p 3 D 2 ), so for the relevant parameters the number of or¬ 
bitals n and the bond dimension D the additional cost scale 
as 0(D 3 + n 3 D 2 ) which is lower in cost than the two-site 
DMRG step with 0(n 2 D 3 + n 3 D 2 ). 


Two-site correlations in the different physical bases 

In order to visualize the different physical bases we present 
below plots for the two site correlations present in the ground 
state approximations within the corresponding basis. The two 
site correlation is measured here by the mutual information 
I(q,r) = S 1 ^) + 5' 1 ( r ) ~ S ,1 ( 9 ,r) with q,r £ [n] and 
5 ll (/) denotes the von Neumann entropy of the reduced state 
pi = tr[ n i\ / |r/))(^|. We present the correlation patterns of the 


final states of the individual calculations shown in Fig. 2 in 
the main text. Fig. 00 show the mutual information of the 
final states obtained at D max = 256 in the initial and opti¬ 
mised basis starting in the HF and localised orbitals. Fig. [8][9] 
display the mutual information patterns obtained in the calcu¬ 
lation with bounded truncation error e t rc = 10~ 6 . 

The correlation patterns obtained in the localised (Fig. [6]) 
and optimised bases (Fig. 000 show very similar features, 
highlighting the fact that the localised basis in this system will 
be close to optimal. Note, however, that although the correla¬ 
tions for the localised and optimised version of the localised 
basis are almost the same (compare Fig.[6]and Fig. [7j, the lat¬ 
ter allows for a more efficient approximation of the ground 
state using a QC-DMRG; displaying the fact that the mutual 
information does not contain the full information about the 
optimality of the basis. 


”i H 

■ II 

1 ■ BB 
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cncvj^ocoi-r^oirjpO'^-coT-Lni-h'COCMcocoojo^-co 
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FIG. 4. Mutual information present in the ground state approxima¬ 
tion within the HF basis by an 17(1) x I/(l)-symmetric MPS with 
Dmax = 256 (dark blue dashed curve in the right panel of Fig. 2 in 
the main text). 
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FIG. 5. Mutual information present in the ground state approxima¬ 
tion within the optimised HF basis by an U(l) x I/(l)-symmetric 
MPS with D max = 256 (light blue curve in the right panel of Fig. 2 
in the main text). 
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FIG. 6. Mutual information present in the ground state approxima¬ 
tion within the localised basis by an (7(1) x ?7(l)-symmetric MPS 
with D mm = 256 (red dashed-dotted curve in the right panel of Fig. 
2 in the main text). 
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FIG. 7. Mutual information present in the ground state approxi¬ 
mation within the optimised localised basis by an (7(1) x (7(1)- 
symmetric MPS with D mix = 256 (orange curve in the right panel 
of Fig. 2 in the main text). 
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FIG. 8. Mutual information present in the ground state approxima¬ 
tion within the HF basis by an (7(1) x U (l)-symmetric MPS with a 
bounded truncation error e t rc < 10~ 6 (dark blue dashed curve in the 
left panel of Fig. 2 in the main text). 
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tion after the tenth basis optimisation by an 17(1) x U (l)-symmetric 
MPS with a bounded truncation error e t rc < 1CP 6 (light blue curve 
in the left panel of Fig. 2 in the main text). 

















































































